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1. Introduction 

The purpose of this note is to illustrate the program MINT, a FORTRAN program for Monte 
Carlo adaptive integration and generation of unweighted distributions. The program per- 
forms the following task: given a function f(x 1 , . . . , x n ), defined in the unit n-dimensional 
cube, and such that the function 

f{x\...,x p ) = f 1 dx p+1 f 1 dx p+2 ... f 1 dx n f(x\...,x n ) (1.1) 
Jo Jo Jo 

is positive for x 1 , . . . , x p in the p-dimensional unit cube, it computes the integral of /, and 
generates x 1 , . . . , x p points in the unit cube, distributed with probability 

/(x 1 ,...,aJ')dx 1 ...dx p . (1.2) 

The case p = n is contemplated, i.e. MINT can also perform the task of unweighting a 
positive distribution. 

A popular program to perform adaptive Monte Carlo multi-dimensional integration is 
the VEGAS program [1]. It uses an importance-sampling method, where the sampling rate 
is independent for each coordinate. The method works well for factorizable singularities, 
and it has the advantage that the sampling information one needs to store grows only 
linearly with the number of dimensions. 

A common problem encountered in particle physics phenomenology, besides the in- 
tegration of a given multidimensional function, is the generation of its arguments with a 
probability proportional to its value. The SPRING-BASES program [2] use the VEGAS 
algorithm to perform the integration of a positive function, and then can generate events 
distributed with a probability proportional to the integrand. It first stores the integration 
result and the maximum value of the function for each cell of the adaptive mesh found by 
VEGAS. In the generation stage, a cell is chosen with a probability proportional to the 
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corresponding value of the integral, and then a point in the cell is generated using the hit 
and miss technique. This method is highly efficient, but it has the disadvantage that the 
amount of storage space it requires grows exponentially with the dimension. 

MINT is a replacement for the SPRING-BASES package. It differs from it in the 
following way. In order to generate the phase-space point, it does not store the value of 
the integral for each cell. It only stores an upper bound of the value of the function in each 
cell. The multidimensional stepwise function that equals the upper bound of the function 
to be integrated in each cell is in fact an upper bound for the whole function, and it is easy 
to generate phase space points distributed according to it. Using again the hit and miss 
technique, we can then generate points according to the original distribution. In order 
to save space, MINT uses an upper bounding functions that is the products of a set of 
step-wise functions, each of them associated with a coordinate. The storage requirement 
for such a function grows only linearly with the dimension. The ability to deal with non- 
positive functions is achieved in MINT by folding over (in a sense that will be specified later) 
the integrand along the directions p + 1, . . . , n. 

The problem of unweighting a distribution of the form (1.1) arises in the context of 
the method of refs. [3, 4, 5] for the inclusion of next-to-leading corrections in Monte Carlo 
generated events. 



2. The algorithm 

The algorithm is in essence the VEGAS algorithm. It is implemented in MINT in the follow- 
ing way. Assume that we deal with the n-dimensional integral of a function /(x 1 , . . . , x n ) ^ 
in the unit hypercube. We divide the [0, 1] interval for each coordinate in m bins of vari- 
able length, 

Xi_i ^ x k ^ xf for I = 1, . . . m and k = 1, . . . , n. 

We then define n monotonic, continuous functions h k (y k ), with < y k < 1 and k = 1, . . . , n, 
such that 

,fc I M k 



m 



h K — = x f for/ = 0,...,m, (2.1) 



and linear (i.e. having constant first derivative) in all the intervals (I — l)/m < y < l/m. 
We have 



/ f(x)d n x = J f {h{y)) f[ (2.2) 



k=i 

,k 



Observe that, given y , if I is the bin where y lies (i.e. (Z — l)/m < y < l/m), we have 

dh k (y k ) 



dy k 



= {xf - x k _i) x m. (2.3) 



We begin the adaptation process with xf = l/m for all k = 1, . . . , n, I = 1, . . . , m. We 
would like to find optimal h k functions, such that the integration process has small errors. 
This is done as follows. We perform several iteration of the integration. In each iteration 
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we generate a set Xjq of N random points y in the unit hypercube O^y^l, k = 1, . . . .,n, 
and estimate the integral with the formula 

'-vE /<%» n^- <«> 

fe=l 

To each y we associate a point x = h(y), and we call l k the bin where x k lies. We accumulate 
the value of / and the number of hits in two arrays R k and N k , with k = 1, . . . ,n and 
I = 1, . . . , m 

n 

R k = J2e(i-i^<i)x f{x) J] (4' - xfLj, jv* = E^-i^y fe < 0, (2-5) 

j/gXjv fe'=i j/eXjv 

where 

n 

fc=l 

is the volume of the cell where x lies. After all the N points have been generated, we define 
for all k = 1, . . . ., n and I = 0, . . . , m 

^E|^ = ° ( 2 - 6 ) 
j=i j 

and a corresponding continuous piecewise linear function i k {x k ) such that i k {x k ) = for 

i = 0, , m. Notice that i k {x k ) is an estimate of the integral of / in the n — 1 dimensional 

hyper-rectangle defined by the condition ^ x k < x k . We now find the new x k points by 
solving the equation 

i*(x*) = i fc (l)- (2.7) 
m 

and then we go to the next iteration. It is clear that when all the Rj/N k are equal, the 
procedure has reached stability, the xf no longer change, and we are probing regions of 
equal importance with (in the average) the same number of points. Furthermore, this 
procedure is nearly optimal if 

f(x\....,x n )=f 1 (x 1 )x....xf n (x n ). (2.8) 
3. Folded integration 

If the given function / is not positive definite, but its integral over a subset of the integration 
variables is positive, we can turn it into a positive function by folding it over itself a sufficient 
number of times in the given subset of the integration variables. More precisely, we do the 
following. Suppose we want to fold each k coordinate times. We define the function 



-i c -i f it 



Pi 1 Pn 

tl=U l n =0 



f (Hy)) n 



dh k {y K 



dy k 
k=i y 



(3.1) 

y=y{z,i) 
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where 



y k (z k ,l k )-- " + ' 



Pk 

so that, as the variable z k spans the [0, 1] interval, y k spans each of thepfc equal subintervals 

^ < y k < (3.2) 
Pk Pk 

It is clear that 

J f(x)d n x = J f( h (y))p^^^ldy k = J f(z\...,z n )dz 1 ...dz n . (3.3) 

It is convenient for practical purposes to choose the p^ among the integer divisors of m. 
4. Generation of p-tuples 

In order to generate p-tuples x 1 , . . . , x p , with p ^ n, and a probability distribution 

f f(x)dxP +1 ...dx n , 

1 r » x , n forp < n 

p(x\...,xn = l f { x f i x)dx , , (4.i) 

^Jifek forp = fi 

we generate the whole n-tuple, and then discard the p + 1, . . . , n coordinates. We can use 
an arbitrary amount of folding on the p + 1, . . . , n coordinates, while we must keep the 
1, . . . ,p coordinates unfolded (i.e. p 1 , . . . ,p q = 1). We then proceed with the generation of 
the z variables using the / function, with probability 

We seek for / an upper bound of the form 

f(z 1 ,...,z n )^u 1 (z 1 )x...xu n (z n ), (4.3) 

where u k {z k ) are stepwise functions in the < z k < 1 interval divided into m/pk equal 
subintervals. The u k {z k ) are initialized according to 



i l/n 

u k (z k ) 



J f(z\...,z n ) dz 1 ...dz r - 



(4.4) 



We perform a large number of calls to the function /, at (uniform) random values of its 
argument. If the bound is violated 

f(z\...,z n )^u 1 (z 1 ) x ... xu n (z n ), (4.5) 

each of the u k (zk) is increased by a fixed factor / in the subinterval containing z^. The 
value 

/-i + iiS («) 

is found to work reasonably well. After a sufficiently large number of calls, the u will 
stabilize. 

In practice, in the MINT program, the upper bounding envelope is computed during 
the folded integration. 
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5. The code 



The code is available at the URL 

http : //moby .mib . inf n. it/~nason/POWHEG/FNOpaper/mint-integrator . f . 

It is composed by two user-callable routines, mint and gen. Furthermore, the common 
block 

integer if old(ndimmax) 
common/ cif old/ if old 

must be set, specifying how many times each dimension is folded. The user function must 
have the form 
function fun(x,w,ifl) 
real * 8 fun,w,x(ndim) 

When called with if 1=0 it must return fun=f (x)*w. When called with if 1=1 it returns 
fun=f (x)*w, but it may assume that, since the previous if 1=0 call the variables that are 
not folded, and all the values that have been computed with them, have remain unchanged. 
The return value f (x)*w is not used by mint in this case. When called with if 1=2, fun 
must return the sum of all the return values up to (and including) the last if 1=0 call. This 
apparently cumbersome procedure is needed to allow for enough flexibility, in order for a 
program to be able to compute the fraction of positive and negative contributions to the 
folded integral, which is needed in POWHEG applications. 
To run the program, one first calls the subroutine 

real * 8 xgrid(50,ndimmax) ,xint,ymax(50,ndimmax) ,ans,err 
integer ndim , ncallsO , nitmax , imode 

call mint (f un , ndim , ncallsO , nitmax , imode , xgr id , xint , ymax , ans , err) 
c ndim: dimension of x in fun (ndim<=ndimmax) 
c ncallsO: maximum number of calls per iteration 
c nitmax: number of iterations 
c imode : flag 

• When called with imode=0, mint performs the integration of the absolute value of 
the function, finds the optimal grid xgr id, stores the answer in xint and ans, and 
the error in err. 

• When called with imode=l, mint performs the folded integration. The grid is kept 
fixed at this stage. The array 

integer if old(ndimmax) 
common/ cif old/ if old 

controls the folding, i.e. ifold(k) is the number of folds for x(k). The number of 
folds must be an integer divisor of the number of bins, which is fixed to 50. The array 
xgrid and the value xint must have already been filled by a previous call to mint 
with if 1=0. The upper bounding envelope of the folded function is also computed 
at this stage, and stored in the array ymax. 
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The value and error of the integral are returned in ans and error. Once ymax has been 
setup in this way, one can call gen to generate events with the distribution of the folded 
function. Notice that negative values of the called function are not allowed at this stage. 
One calls 
imode=0 

call gen (fun , ndim , xgrid , ymax , imode , x) 

imode=l 

do j=l, 10000 

call gen (fun , ndim , xgrid , ymax , imode , x) 

enddo 
imode=3 

call gen (fun , ndim , xgrid , ymax , imode , x) 

where the call with imode=0 initializes the generation, imode=l generate the ndim-tuples, 
and imode=3 prints out generation statistics. After a call to gen with imode=l, one can 
assume that the last call to the function fun was performed with the generated value of 
x, so that parameters depending upon x that are stored in common blocks by fun have 
consistent values. 
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